from sympy import *
n , k = symbols('n, k')
c6 = Function('c6')
c5 = Function('c5')
f0 = lambda n, k: 1
f1 = lambda n, k: k
f2 = lambda n, k: k * f1(n, k) - 1 * n * f0(n, k)
expand(f2(n, k))
f3 = lambda n, k: k * f2(n, k) - 2 * n * f1(n, k)
expand(f3(n, k))
f4 = lambda n, k: k * f3(n, k) - 3 * n * f2(n, k) + 2 * n
expand(f4(n, k))
expand((f4(n+1, k+1) + f4(n+1, k-1))/2 )
f5 = lambda n, k: k * f4(n, k) - 4 * n * f3(n, k) + 8 * n * k
expand(f5(n, k))
expand((f5(n+1, k+1) + f5(n+1, k-1))/2 )
f6 = lambda n, k: k * f5(n, k) - 5 * n * f4(n, k) + 20 * n * k ** 2 - 16 * n - 20 * n ** 2
expand(f6(n, k))
expand((f6(n+1, k+1) + f6(n+1, k-1))/2 ) - expand(f6(n, k))
simplify(20 * n * k ** 2 - 36 * n - 20 * n ** 2 + 20 * n)
2 * n * f0(n,k)
8 * n * f1(n, k)
simplify(20 * n *(f2(n, k)) - 16 * n * f0(n,k))
f7 = lambda n, k: k * f6(n, k) - 6 * n * f5(n, k) + 40 * n * f3(n, k) - 96 * n * f1(n, k)
expand(f6(n, k))
expand((f7(n+1, k+1) + f7(n+1, k-1))/2 ) - expand(f7(n, k))
simplify(40 * n * f3(n, k) - 96 * n * f1(n, k))
f8 = lambda n, k: k * f7(n, k) - 7 * n * f6(n, k) + 70 * n * f4(n, k) - 336 * n * f2(n, k) + 272 * n * f0(n, k)
expand((f8(n+1, k+1) + f8(n+1, k-1))/2 ) - expand(f8(n, k))
f9 = lambda n, k: k * f8(n, k) - 8 * n * f7(n, k) + 112 * n * f5(n, k) - 896 * n * f3(n, k) + 2176 * n * f1(n, k)
expand((f9(n+1, k+1) + f9(n+1, k-1))/2 ) - expand(f9(n, k))
f10 = lambda n, k: k * f9(n, k) - 9 * n * f8(n, k) + 168 * n * f6(n, k) - 2016 * n * f4(n, k) + 9792 * n * f2(n, k) - 7936 * n * f0(n, k)
expand((f10(n+1, k+1) + f10(n+1, k-1))/2 ) - expand(f10(n, k))
f11 = lambda n, k: k * f10(n, k) - 10 * n * f9(n, k) + 240 * n * f7(n, k) - 4032 * n * f5(n, k) + 32640 * n * f3(n, k) - 79360 * n * f1(n, k)
expand((f11(n+1, k+1) + f11(n+1, k-1))/2 ) - expand(f11(n, k))
import functools
def f(m):
if m < 0:
return lambda n, k: 0
if m == 0:
return lambda n, k: 1
elif m == 1:
return lambda n, k: k
else:
def alternating_tangent(n):
# https://mathworld.wolfram.com/TangentNumber.html
return 2 ** (2*n) * (2**(2*n)-1) * Abs(bernoulli(2*n)) / (2*n)
def R(n, k):
running_sum = 0
for j in range(1, round(m/2) + 1):
i = m - 2*j
running_sum += (-1)**(j) * alternating_tangent(j) * binomial(m-1, (2 * j - 1)) * f(i)(n, k)
return running_sum
def fm(n, k):
return k * f(m-1)(n,k) + n * R(n, k)
return fm
print(latex(expand(f(8)(n,k))))
k^{7} - 21 k^{5} n + 105 k^{3} n^{2} + 70 k^{3} n - 105 k n^{3} - 210 k n^{2} - 112 k n
f(6)(6,6)
expand(f(4)(n+1,k))
expand(f(4)(n,k) - 6 * f(2)(n, k) + 5*f(0)(n, k))
expand(f(5)(n+1,k))
expand(f(5)(n,k) -10 * f(3)(n,k) + 25 * f(1)(n,k))
expand(f(6)(n+1,k))
expand(f(6)(n,k) - 15* f(4)(n,k) + 75 * f(2)(n, k) - 61 * f(0)(n,k))
# the above coefficients look like http://www.luschny.de/math/seq/SwissKnifePolynomials.html
expand(f(4)(n,k-1) - f(4)(n,k+1))
expand(f(4)(n,k) - 4 * f(3)(n, k) + 6*f(2)(n, k) - 4*f(1)(n, k) + 1 * f(0)(n, k))